function priceapprox=constprice
global x r T t rho K % price parameters
global az bz mz % Zt, OU process parameters
global ay by my % Yt, OU process parameters
global eps v% = 1/ay
global sigma1sqhd sigma1sqbar
global Npath betaK

x = 1; r = 0.1; T = 0.5; t = 0; rho = 0.5; K=0; % x has to be nonzero in order P1 != 0 
az = 1;   bz = 1; mz=1; % let Zt be the simplest OU proess
ay = 100; by = 1; my=1; % Yt 's coeff as OU process, ay = 1/eps;
eps = 1/ay; v = sqrt(by^2/2/ay); % page 68
sigma1sqhd = @(s) exp(s*2); % sigma1 = s
sigma1sqbar = barsigmasq; % <sigma1sq> = <s^2> = v^2+my^2 see page 68 page 86

sigma1sqbar = 20;

Npath = 1; betaK=20; 

sampleprice0=zeros(Npath,1); %sampleprice1=zeros(Npath,1);
%[v2,v3] = coeffv3;

for ii=1:Npath
    [P0,ptptP0,ptptptP0] = conditionedprice;
    sampleprice0(ii) = P0;
%    v2bar = -v2*intZt(ii); v3bar = -v3*intZt(ii);
%    sampleprice1(ii) = v2bar*x^2*ptptP0+v3bar*x^3*ptptptP0;
    disp(ii)
end

price0 = sum(sampleprice0)/Npath;
%price1 = sum(sampleprice1)/Npath;

%priceapprox  = price0+price1;
priceapprox = price0;
end


function [P0,ptptP0,ptptptP0] = conditionedprice
global x r T t K sigma1sqbar

mu = r*(T-t)-sigma1sqbar/2*(T-t);
sigma = sqrt(sigma1sqbar*(T-t));

P0 = quadgk(@(s) s.*lognpdf(s/x,mu,sigma)/x,-Inf, Inf);
ptptP0 = quadgk(@(s) ((mu-log(s)).^2/sigma^4-1/sigma^2).*max(s-K,0).*lognpdf(s,mu,sigma),-Inf, Inf);
ptptptP0 = quadgk(@(s) ((mu-log(s)).^3/sigma^6+3*(mu-log(s))/sigma^4).*max(s-K,0).*lognpdf(s,mu,sigma),-Inf, Inf);

%P0 = exp(-r*(T-t))*P0;
ptptP0 = exp(-r*(T-t))*ptptP0;
ptptptP0 = exp(-r*(T-t))*ptptptP0;
end

function [v2,v3] = coeffv3
global rho ay v

phiprime = @newapproach;
average=quadgk(@(y) phiprime(y).*y, -Inf, Inf); 
% here phi' is actually phi'(y)*normpdf(y,m,v), so average actually
% = y*phi'(y)*normpdf(y,m,v)=<yphi'(y)>

v3 = v/sqrt(2*ay)*rho*average;
v2 = 2*v3;
end

function result=newapproach(y) % the verification see ../testing/testcoeff3.m
global v my sigma1sqbar
z=y-my;
%A1=@(z) (z<0).*(sqrt(pi)/4*erfc(-z)-z/2.*exp(-z.^2))+(z>=0).*(sqrt(pi)/2-z.*exp(-z.^2)/2-sqrt(pi)/4*erfc(z));
A1=@(z) sqrt(pi)/2-z.*exp(-z.^2)/2-sqrt(pi)/4*erfc(z);
result = 1/sqrt(2*pi)/v*(sqrt(2*v^2)^3*A1(z/sqrt(2)/v)-2*v^2*my*exp(-z.^2/2/v^2)+my^2*sqrt(2*pi*v^2)/2*erfc(-z/sqrt(2*v^2)))-sigma1sqbar/2*erfc(-z/sqrt(2*v^2));
result = result/v^2;
% result <==> quadgk(@(s) (sigma1(s)-sigma1sqbar).*normpdf(s,my,v), -Inf, y)      /v^2;
% which is phi'(y)*normpdf(y,m,v);

end

function y = barsigmasq
global my v sigma1sqhd
%y = quadgk(@(s) sigma1sqhd(s).*normpdf(s,my,v), -Inf, Inf);
y = quadgk(@(s) 1/sqrt(2*pi)/v*exp(2*s-(s-my).^2/2/v^2), -Inf, Inf);
%y = v^2+my^2;
end